#include "../PhysicalField/PhysicalField.h"

#include <limits>
#include <iostream>

#include "../Input/PicParams.h"
#include "../Species/Species.h"
#include "../Field/Field.h"
#include "../FieldSolver/SolverFactory.h"

using namespace std;


//---------------------------------------------------------------------------------------------------------------------
//Constructor for the virtual class PhysicalField
//---------------------------------------------------------------------------------------------------------------------
PhysicalField::PhysicalField(PicParams* params_, ParallelMPI* pmpi_) :
    params(params_),
    pmpi(pmpi_),
    n_species(params_->species_param.size())
{
    species_param = params->species_param;

    dims_total = 1;
    dims_totalglobal = 1;
    for (int i = 0 ; i < params->n_dim_field; i++)
    {
        dims_total    *= params->dims_field[i];
        dims_totalglobal *= params->dims_field_global[i];
    }

    Ex_=NULL;
    Ey_=NULL;
    Ez_=NULL;
    Bx_=NULL;
    By_=NULL;
    Bz_=NULL;
    rho_=NULL;

    Ex_avg=NULL;
    Ey_avg=NULL;
    Ez_avg=NULL;
    Bx_avg=NULL;
    By_avg=NULL;
    Bz_avg=NULL;

    rho_s.resize(n_species);
    rho_s_avg.resize(n_species);
    rho_s_global.resize(n_species);
    rho_s_global_avg.resize(n_species);

    Vx_s.resize(n_species);
    Vx_s_avg.resize(n_species);
    Vx_s_global.resize(n_species);
    Vx_s_global_avg.resize(n_species);

    Vy_s.resize(n_species);
    Vy_s_avg.resize(n_species);
    Vy_s_global.resize(n_species);
    Vy_s_global_avg.resize(n_species);

    Vz_s.resize(n_species);
    Vz_s_avg.resize(n_species);
    Vz_s_global.resize(n_species);
    Vz_s_global_avg.resize(n_species);

    Vp_s.resize(n_species);
    Vp_s_avg.resize(n_species);
    Vp_s_global.resize(n_species);
    Vp_s_global_avg.resize(n_species);

    T_s.resize(n_species);
    T_s_avg.resize(n_species);
    T_s_global.resize(n_species);
    T_s_global_avg.resize(n_species);

    Tx_s.resize(n_species);
    Tx_s_avg.resize(n_species);
    Tx_s_global.resize(n_species);
    Tx_s_global_avg.resize(n_species);

    Ty_s.resize(n_species);
    Ty_s_avg.resize(n_species);
    Ty_s_global.resize(n_species);
    Ty_s_global_avg.resize(n_species);

    Tz_s.resize(n_species);
    Tz_s_avg.resize(n_species);
    Tz_s_global.resize(n_species);
    Tz_s_global_avg.resize(n_species);

    q_s.resize(n_species);
    q_s_avg.resize(n_species);
    q_s_global.resize(n_species);
    q_s_global_avg.resize(n_species);

    for(int ispec=0; ispec<n_species; ispec++)
    {
        rho_s[ispec]            = NULL;
        rho_s_avg[ispec]        = NULL;
        rho_s_global[ispec]     = NULL;
        rho_s_global_avg[ispec] = NULL;

        Vx_s[ispec]             = NULL;
        Vx_s_avg[ispec]         = NULL;
        Vx_s_global[ispec]      = NULL;
        Vx_s_global_avg[ispec]  = NULL;

        Vy_s[ispec]             = NULL;
        Vy_s_avg[ispec]         = NULL;
        Vy_s_global[ispec]      = NULL;
        Vy_s_global_avg[ispec]  = NULL;

        Vz_s[ispec]             = NULL;
        Vz_s_avg[ispec]         = NULL;
        Vz_s_global[ispec]      = NULL;
        Vz_s_global_avg[ispec]  = NULL;

        Vp_s[ispec]             = NULL;
        Vp_s_avg[ispec]         = NULL;
        Vp_s_global[ispec]      = NULL;
        Vp_s_global_avg[ispec]  = NULL;

        T_s[ispec]              = NULL;
        T_s_avg[ispec]          = NULL;
        T_s_global[ispec]       = NULL;
        T_s_global_avg[ispec]   = NULL;

        Tx_s[ispec]              = NULL;
        Tx_s_avg[ispec]          = NULL;
        Tx_s_global[ispec]       = NULL;
        Tx_s_global_avg[ispec]   = NULL;

        Ty_s[ispec]              = NULL;
        Ty_s_avg[ispec]          = NULL;
        Ty_s_global[ispec]       = NULL;
        Ty_s_global_avg[ispec]   = NULL;

        Tz_s[ispec]              = NULL;
        Tz_s_avg[ispec]          = NULL;
        Tz_s_global[ispec]       = NULL;
        Tz_s_global_avg[ispec]   = NULL;

        q_s[ispec]              = NULL;
        q_s_avg[ispec]          = NULL;
        q_s_global[ispec]       = NULL;
        q_s_global_avg[ispec]   = NULL;
    }

}



PhysicalField::~PhysicalField()
{
    delete Ex_;
    delete Ey_;
    delete Ez_;
    delete Bx_;
    delete By_;
    delete Bz_;
    delete rho_;

    if (Ex_avg!=NULL)
    {
        delete Ex_avg;
        delete Ey_avg;
        delete Ez_avg;
        delete Bx_avg;
        delete By_avg;
        delete Bz_avg;
    }

    for (int ispec=0; ispec<n_species; ispec++)
    {
        delete rho_s[ispec];
    }


}


void PhysicalField::initRhoJ(vector<Species*>& vecSpecies, Projector* Proj)
{


}



void PhysicalField::incrementAvgFields(int itime)
{
    if( (itime - 1) % params->timestep_io == 0 )
    {
        rho_global_avg->set_value(0.0);
        phi_global_avg->set_value(0.0);
        Ex_global_avg->set_value(0.0);

        for(int ispec=0; ispec<n_species; ispec++)
        {
            rho_s_avg[ispec]->set_value(0.0);
        }
    }

    //Calculate the sum values for global rho phi Ex and Ey
    if( (itime % params->timestep_io) > (params->timestep_io - params->timestep_average) || (itime % params->timestep_io) == 0 )
    {
        for(int i = 0 ; i < dims_totalglobal ; i++)
        {
            (*rho_global_avg)(i) += (*rho_global)(i);
            (*phi_global_avg)(i) += (*phi_global)(i);
            (*Ex_global_avg)(i)  += (*Ex_global)(i);
        }

        //Calculate the sum values for density of each species
        for(int ispec = 0; ispec < n_species; ispec++)
        {
            //all fields are defined on the primal grid
            for(int i = 0 ; i <dims_total; i++)
            {
                (*rho_s_avg[ispec])(i) += (*rho_s[ispec])(i);
            }
        }
    }

    //calculate the averaged values
    if( itime % params->timestep_io == 0 )
    {
        for(int i=0 ; i<dims_totalglobal ; i++)
        {
            (*rho_global_avg)(i) /= params->timestep_average;
            (*phi_global_avg)(i) /= params->timestep_average;
            (*Ex_global_avg)(i)  /= params->timestep_average;
        }
        for(int ispec = 0; ispec < n_species; ispec++)
        {
            for(int i = 0 ; i < dims_total ; i++)
            {
                (*rho_s_avg[ispec])(i) /= params->timestep_average;
            }
        }
    }

}


void PhysicalField::restartRhoJ()
{
    rho_->set_value(0.0);
}

void PhysicalField::restartRhoJs(int ispec)
{
    rho_s[ispec]->set_value(0.0);
}


void PhysicalField::computeTotalRhoJ()
{
    for(int ispec=0; ispec<n_species; ispec++)
    {
        for(int i = 0 ; i < dims_total; i++)
        {
            (*rho_)(i) += ( (*rho_s[ispec])(i) * species_param[ispec].charge );
        }
    }
}


void PhysicalField::gatherFields(ParallelMPI *pmpi)
{
    for(int i = 0; i < rho_s.size(); i++)
    {
        pmpi->gatherField1(rho_s_global[i], rho_s[i]);
    }
    pmpi->gatherField1(rho_global, rho_);
}


void PhysicalField::gatherAvgFields(ParallelMPI *pmpi)
{
    for(int i = 0; i < rho_s.size(); i++)
    {
        pmpi->gatherField1( rho_s_global_avg[i], rho_s_avg[i] );
    }
}
